Packages

library(plyr)
library(tidyverse)
library(popdemo)
library(popbio)

Load functions.

source('R/pop01_param_poun.R')
source('R/pop03_doproj.R')
source('R/pop03_doproj_stoch.R')

Projections

Population projection parameters.

Each row holds a unique set of parameters that reflect diverse recovery (conservation action) and extinction (threat) scenarios.

dt <- pop01_param_poun()
 mydigits <- c(NA, NA, 1, 0, 0, 0, 2, 2, 2, 0, 0, 0, 2, 2, 0, 0 ,0, 2, 2, NA)
slice_head(dt, n = 5)  |> 
  kableExtra::kbl(digits = mydigits) |>
  kableExtra::kable_styling(full_width = F,  latex_options = "hold_position")
species type first_year a1 a2 a3 a4 b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 d3 d4 akey
Podocnemis unifilis base 0.0 0 0 0 10.74 0.0 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0
Podocnemis unifilis base 0.1 0 0 0 10.74 0.1 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.1
Podocnemis unifilis base 0.2 0 0 0 10.74 0.2 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.2
Podocnemis unifilis base 0.3 0 0 0 10.74 0.3 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.3
Podocnemis unifilis base 0.4 0 0 0 10.74 0.4 0.33 0 0 0 0.17 0.29 0 0 0 0.11 0.93 poun_base_0.4

Projection models.

Deterministic models

nf <- 100
dt$adultF_n <- nf
# project
dout <- plyr::ddply(dt, 
                    c("species", "type", "first_year","akey"), .fun =  pop03_doproj)
dout$arun <- 1

# Model summaries
model_sum <- dout |> 
  group_by(species, type, first_year, lambda, gen_time) |> 
  summarise(fem_t0 = max(fem_t0), 
            fem_min = min(fem),
              fem_max = max(fem)) |> 
  ungroup()
lambda_n <- length(unique(model_sum$lambda)) # 50
lambda_mean <- mean(model_sum$lambda) # 0.9432
lambda_sd <- sd(model_sum$lambda) # 0.1506
lambda_min <- min(model_sum$lambda) # 0.4659
lambda_max <- max(model_sum$lambda) # 1.1539

Stochastic models

# Stochastic
#data frame with runs for processing
#nruns <- 100 # 100 gives same pattern as 50
nruns <- 50
dt_stoch <- dt[rep(seq_len(nrow(dt)), nruns), ]
dt_stoch$arun <- rep(1:nruns, each = 50)
# Approx 7 minutes. 1 million rows.
dout_stoch <- plyr::ddply(dt_stoch, 
                    c("arun", "species", "type", "first_year","akey"), 
                    .fun =  pop03_doproj_stoch)
table(dout_stoch$model)

Combine results.

# Combine data for plotting
dout_all <- dplyr::bind_rows(dout |> dplyr::select(arun, model, type, first_year, 
                                                   akey, ayear, lambda,
                                            fem, fem_t0, fem_diff, change50_flag, 
                                            double_flag), 
                             dout_stoch |> dplyr::select(arun, model, type, first_year, 
                                                         akey, ayear, 
                                                         lambda,
                                            fem, fem_t0, fem_diff, change50_flag, 
                                            double_flag))
# Limit adult female number to maximum (10 x original).
dout_all[which(dout_all$fem > (nf * 10)), 'fem' ] <- (nf * 10)
# Factors in right order
dout_all$modelf <- 1
dout_all[which(dout_all$model=="Stochastic uniform") , 'modelf'] <- 2
dout_all[which(dout_all$model=="Stochastic equal") , 'modelf'] <- 3
dout_all[which(dout_all$model=="Stochastic bad x2") , 'modelf'] <- 4
dout_all[which(dout_all$model=="Stochastic bad x4") , 'modelf'] <- 5
dout_all$modelf <- as.factor(dout_all$modelf)
levels(dout_all$modelf) <- c("Deterministic", "Stochastic uniform", 
                          "Stochastic equal", "Stochastic bad x2", 
                          "Stochastic worst x4")
unique(dout_all$modelf)
table(dout_all$modelf)
dout_all$typef <- 1
dout_all[which(dout_all$type=="female-hunt 2.5%") , 'typef'] <- 2
dout_all[which(dout_all$type=="female-hunt 10%") , 'typef'] <- 3
dout_all[which(dout_all$type=="female-hunt 25%") , 'typef'] <- 4
dout_all[which(dout_all$type=="female-hunt 50%") , 'typef'] <- 5
dout_all$typef <- as.factor(dout_all$typef)
levels(dout_all$typef) <- c("base", "female-hunt 2.5%", 
                          "female-hunt 10%", "female-hunt 25%", 
                          "female-hunt 50%")
table(dout_all$typef)
# first year survival
dout_all$first_yearf <- as.factor(dout_all$first_year)
fylev <- paste("first-year\nsurvival\n", seq(0, 0.9, by = 0.1), sep = "")
levels(dout_all$first_yearf) <- fylev
# Export for future use
saveRDS(dout_all, "inst/other/dout_all.rds")

Model summaries.

# load data
dout_all <- readRDS("inst/other/dout_all.rds")

# Make unique model ID. boot mean is same as mean (at least to 6 decimal places)
# 250 projection models.
model_ref <- dout_all |>
  group_by(akey, modelf, typef, first_year, first_yearf, arun) |> 
  summarise(lambda = min(lambda), 
            yc = length(unique(ayear))) |> 
  ungroup() |> 
  group_by(akey, modelf, typef, first_year, first_yearf) |>
  summarise(count_runs = length(unique(arun)), 
            count_years = min(yc),
            lambda_mean = mean(lambda), 
            lambda_min = min(lambda), 
            lambda_max = max(lambda), 
            lambda_boot_lcl = Hmisc::smean.cl.boot(lambda)["Lower"],
            lambda_boot_ucl = Hmisc::smean.cl.boot(lambda)["Upper"]
            ) |> 
  ungroup() |> 
  arrange(typef, first_year, modelf) |>
  mutate(modelid = paste(akey, as.numeric(modelf), sep = "_"), 
         modelkey = row_number()) |> 
  relocate(modelkey, modelid)
## `summarise()` has grouped output by 'akey', 'modelf', 'typef', 'first_year',
## 'first_yearf'. You can override using the `.groups` argument.
## `summarise()` has grouped output by 'akey', 'modelf', 'typef', 'first_year'.
## You can override using the `.groups` argument.
# summaries of  10050 runs.
#model_all_sum <- dout_all |> 
#  group_by(akey, modelf, typef, first_year, arun) |> 
#  summarise(count_years = length(ayear), 
#            count_years_unique = length(unique(ayear))) |> 
#  ungroup()

Plots

mean_values <- model_ref |> 
  group_by(typef) |> 
  summarise(mean_l = mean(lambda_mean))


fig_model_lambda <- model_ref |> 
  ggplot(aes(x = first_year, y = lambda_mean, colour = modelf)) + 
  geom_hline(data = mean_values, aes(yintercept = mean_l), 
             linetype = "dashed") + 
  geom_point() + 
  stat_smooth(se = FALSE) + 
  scale_x_continuous("first year survival and graduation", 
                     breaks = c(0, 0.2, 0.4, 0.6, 0.8)) + 
  scale_color_viridis_d() +
  facet_wrap(~ typef, ncol = 5) + 
  theme(legend.position="top") + 
  guides(colour=guide_legend(nrow=2,byrow=TRUE)) + 
  labs(y = "lambda", colour = "model")
# save
png(file = "inst/other/fig_unifilis_lambda.png", bg = "transparent", 
    type = c("cairo"), 
    width = 8, height = 4, units = "in", res=600)
fig_model_lambda
invisible(dev.off())
Check plot.
Model lambda for *Podocnemis unifilis*

Model lambda for Podocnemis unifilis

Projections over time.

fig_proj <- dout_all |> 
  ggplot(aes(x = ayear, y = fem, colour = modelf)) + 
  geom_point(size=0.1, alpha=0.2) + 
  stat_smooth(se = FALSE) + 
  scale_colour_viridis_d("model") + 
  scale_y_continuous(limits = c(0, 1000)) +
  labs(y="Reproductive females",
       x="Time (years)",
       ) +
  theme_bw() +
  theme(plot.title.position = "plot") + 
  facet_grid(first_yearf ~ typef)

# save
png(file = "inst/other/fig_unifilis_projections.png", bg = "transparent", 
    type = c("cairo"), 
    width = 10, height = 8, units = "in", res=600)
fig_proj
invisible(dev.off())
Check plot.
Model scenarios for *Podocnemis unifilis*

Model scenarios for Podocnemis unifilis

Number of adult females along 80 km2 of river.

# Here use density values from "tracaja_dist_5km_4z_beforeafter.R" to establish "impact"
# Adult population before - after
an_b4 <- ceiling((80 * 1.0035150)) # before = 81 in 80 km2 of river
an_b4_lci <- ceiling((80 * 0.38838535))
an_b4_uci <- ceiling((80 * 2.5928949))
an_aft <- ceiling((80 *  0.1542282)) # after  = 13 in 80 km2 of river
an_aft_lci <- ceiling((80 * 0.04021797))
an_aft_uci <- ceiling((80 * 0.5914359))